FFT (Cooley-Tukey) の式変形
式変形をすると、半分で切って重ねたような数列のフーリエ変換として表せる
$ W_Nは$ 1の原始$ N乗根とする
スタート：$ \hat{A}_i = \sum_{j = 0}^{N - 1} A_j W_N^{ij}
偶数項目について
$ \hat{A}_{2i} = \sum_{j = 0}^{N - 1} A_j W_N^{2ij}
$ W_{N/2}を使って表す
$ = \sum_{j = 0}^{N - 1} A_j W_{N/2}^{ij}
前半と後半に分ける
$ = \sum_{j = 0}^{N/2 - 1} A_j W_{N/2}^{ij} + \sum_{j = N/2}^{N - 1} A_j W_{N/2}^{ij}
$ jの範囲を揃える
$ = \sum_{j = 0}^{N/2 - 1} A_j W_{N/2}^{ij} + \sum_{j = 0}^{N / 2 - 1} A_{j + N/2} W_{N/2}^{i(j + N/2)}
$ W_{N/2}^{i(N/2)} = 1より
$ = \sum_{j = 0}^{N/2 - 1} A_j W_{N/2}^{ij} + \sum_{j = 0}^{N / 2 - 1} A_{j + N/2} W_{N/2}^{ij}
前半と後半をくっつけると離散フーリエ変換の式になった
$ = \sum_{j = 0}^{N/2 - 1} \left(A_j + A_{j + N/2}\right) W_{N/2}^{ij}
奇数項目について
$ \hat{A}_{2i + 1} = \sum_{j = 0}^{N - 1} A_j W_N^{(2i + 1) j}
$ W_N^{(2i+1)j}を分解する
$ = \sum_{j = 0}^{N - 1} A_j W_N^{j} W_N^{2ij}
$ W_{N/2}を使って表す
$ = \sum_{j = 0}^{N - 1} A_j W_N^{j} W_{N/2}^{ij}
前半と後半に分ける
$ = \sum_{j = 0}^{N/2 - 1} A_j W_N^{j} W_{N/2}^{ij} + \sum_{j = N/2}^{N - 1} A_j W_N^{j} W_{N/2}^{ij}
$ jの範囲を揃える
$ = \sum_{j = 0}^{N/2 - 1} A_j W_N^{j} W_{N/2}^{ij} + \sum_{j = 0}^{N/2 - 1} A_{j + N/2} W_N^{j + N/2} W_{N/2}^{i(j + N/2)}
$ W_{N/2}^{i(N/2)} = 1より
$ = \sum_{j = 0}^{N/2 - 1} A_j W_N^{j} W_{N/2}^{ij} + \sum_{j = 0}^{N/2 - 1} A_{j + N/2} W_N^{j + N/2} W_{N/2}^{ij}
前半と後半をくっつけると離散フーリエ変換の式になった
$ = \sum_{j = 0}^{N/2 - 1} (A_j + A_{j + N/2} W_N^{N / 2}) W_N^{j} W_{N/2}^{ij}
これを踏まえると、離散フーリエ変換は以下のプログラムにより可能
ここでは数列は複素数からなり、$ W_N = e^{i\frac{2\pi}{N}}, W_N^{N/2} = -1とする
code:rust
use num_complex::*;
use num_traits::*;
fft_internal(A, A.len(), 1, 0, 0);
}
for x in A.iter_mut() {
x = x.conj();
}
fft_internal(A, A.len(), 1, 0, 0);
for x in A.iter_mut() {
x = x.conj() / A.len() as f64;
}
}
fn fft_internal(A: &mut f64, n: usize, s: usize, q: usize, d: usize) { let w_n = Complex::from_polar(1.0, 2.0 * f64::PI / n as f64);
let mut w = Complex::one();
if n > 1 {
for j in 0 .. n / 2 {
w *= w_n;
}
fft_internal(A, n / 2, 2 * s, q, d);
fft_internal(A, n / 2, 2 * s, q + n / 2, d + s);
} else if q > d {
A.swap(q, d);
}
}
参考